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Summary 

A new version of the LEWICE ice accretion computer code has been developed which calculates 
the ice growth on two-dimensional surfaces incorporating the effects of compressibility through 
solution of the Euler equations. The code is modular and contains separate stand-alone program 
elements that create a grid, calculate the flow field parameters, calculate the droplet trajectory 
paths, determine the amount of ice growth, and plot results. This code increases the applicability 
of ice accretion predictions by allowing calculations at higher Mach numbers. The new elements 
of the code are described. Calculated results are compared to experiment for several cases, includ- 
ing a LEWICE example case and a thin airfoil section at a Mach number of 0.58. 

I. Introduction 

Over the past several years the need for ice accretion prediction capabilities has been growing. 
Many aircraft and ice protection system manufacturers have used the NASA developed ice accre- 
tion code, LEWICE 1 . The LEWICE code predicts the growth of ice on 2D surfaces through appli- 
cation of an inviscid panel method module, a particle trajectory calculation module, and a control 
volume energy balance/ice growth calculation module. The use of this code assists the analyst in 
assessment of potential icing hazards due to ice growth on unprotected surfaces and in the appro- 
priate placement of candidate ice protection systems. Several upgrades to this code are currently 
underway to improve the capability of LEWICE. Some of these include; extension to 3D geome- 
tries^, inclusion of a thermal ice protection system modeP, and the addition of a performance deg- 
radation evaluation capability 4 . 

An additional area of potential improvement for the LEWICE code comes from the introduction 
of an alternate flow code calculation method. The current LEWICE code uses the inviscid panel 
method, S24Y, developed by Hess and Smith~\ This approach has the advantage of being a very 
fast code and produces good results over an adequate range of conditions. There are several prob- 
lems associated with the use of this code which have been documented in the LEWICE User’s 
Manual 1 and by other users of the code*\ The S24Y code uses a Karmen-Tsien compressibility 
correction. This helps extend the range of applicability to higher Mach numbers. However, the ba- 
sis for this correction is a small perturbation approximation, hence thick airfoils and high angles 
of attack may not be accounted for using the correction. Additionally, the convective heat transfer 
values, which contribute significantly to the energy balance/ice growth calculation, may be influ- 
enced by this approximation and thus could affect the resulting ice growth for thick airfoils and 
high angle of attack conditions. 



Another area of concern is the evaluation of velocities near the surface at panel edges. The calcu- 
lated velocities at the panel edges can become much larger than actual values, resulting in the pos- 
sibility of unrealistic reversals in water droplet trajectories. This problem has been addressed in 
the LEW1CE code by the creation of an artificial surface located some distance upstream of the 
actual surface, normally 0.2 to 0.6 percent of the chord length. The droplets, now avoiding the 
flow reversals, impinge on this pseudo-surface and the local collection efficiency is calculated 
based on this surface. 

The collection efficiency value is a measure of the amount of water impinging on the surface com- 
pared to the amount of water passing through a plane upstream of the airfoil bounded by the air- 
foil thickness. The local collection efficiency is a measure of how the droplet impingement pattern 
varies as a function of distance along the surface. Thus, since the pseudo- surface approximation 
can change the surface distance measurement in the region near the leading edge, the collection 
efficiency values can differ from those obtained by using the actual surface contours. This has the 
potential for problems when determining the collection efficiency on iced surfaces. 

In order to include more accurate velocity results near the airfoil surface and to include the effects 
of compressibility, the inviscid panel method was replaced by an Euler code. This change also re- 
quires the use of a grid generator. Additionally, the particle trajectory calculation must be altered 
to use interpolated velocities from the grid as opposed to determining the velocities by evaluating 
the contribution of each panel to the potential function at the desired location. The grid generator 
used has the ability to develop body-fitted grids for extreme shapes, such as glaze ice horns. The 
velocity interpolation calculation is incorporated into a subroutine within the particle trajectory 
code module. 

The selection of an Euler code over a full potential code is based on consideration of later incor- 
poration of Navier-Stokes capability, in order to provide performance degradation calculations. 
Euler codes provide two additional capabilities over full potential methods, the ability to accurate- 
ly capture the behavior across shocks and the ability to calculate the development of vorticity in 
the flow. The former is not of interest in ice accretion calculations, however, the latter capability 
can be useful in describing the interaction of vortices developed by the ice accretions with any 
downstream structures of interest. 

The elimination of restrictions relating to small disturbance theory results in the ability to evaluate 
thicker airfoil shapes and higher angles of attack. This is especially useful for sharp-nosed airfoils 
at angle of attack. These airfoils can produce large pressure gradients, resulting in unrealistic val- 
ues of velocity and convective heat transfer if evaluated using the inviscid panel method. The Eul- 
er code has the potential to more accurately model such conditions. 

As a result of the Euler code producing more realistic velocities in the near field, the need for the 
pseudo-surface approximation is eliminated. The particle trajectories are calculated, using the in- 
terpolated velocity values, until they reach the last grid line before the surface. Once the particle 
has crossed this grid line, the droplet path is no longer altered by the flow field and is considered 
to be tangent to the path line calculated in the previous step. The intersection of that path line with 
the body surface is taken as the impingement location. 

Due to the irregular shapes associated with ice formations, multiple stagnation points can develop 
on the iced surface. This can lead to difficulty in selecting the point at which to start the boundary 
layer calculations for the upper and lower surfaces. LEWICE deals with this situation by either 
asking the user to select one of the stagnation points to start the boundary layer calculations or to 
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create another pseudo- surface encompassing the stagnation points and calculating a single new 
stagnation point for that surface. The use of a grid based code can lead to different approaches to 
the solution of this problem. For example, the velocities at grid point locations away from the sur- 
face could be used to determine which grid line divides the upper surface from the lower surface 
and the surface point associated with that grid line could be the starting point for the boundary 
layer calculations. Such an approach could be automated and eliminate the need for user interac- 
tion as is currently the case in LEWICE. Although this capability has not been required in 
LEWICE/E for the current calculations, it is expected that it will be required and thus will be add- 
ed to a future update of the code. 

II. Code Description 

The code consists of four major modules tied together by a simple command-line user interface. 
Introduction of a graphical user interface is planned and should result in greater ease of use. The 
four code modules are a grid generator, an Euler flow solver, a particle trajectory calculation, and 
an energy balance/ice growth calculation. The code modules produce graphs of Cp distributions, 
residual histories, particle paths, collection efficiency distributions, flow field properties at the 
edge of the boundary layer, thermal conditions on the surface and the calculated ice shape. Con- 
tour plots can be produced which show the conditions in the far field. Output listings from each of 
the modules are also produced in order to allow examination of exact numerical values for param- 
eters of interest. Files for geometry and flow solution information are written in PLOT3D format 
in order to allow compatibility with other grid generators and flow solvers which could be substi- 
tuted for those used in this code. 

H-l. Grid Generation Module 

The grid generator used in LEWICE/E is a hyperbolic equation solver developed by Barth . This 
code creates a C-grid around the user defined surface and can be easily applied to complex surfac- 
es such as an iced airfoil. The number of grid points and grid point spacing along the surface can 
be adjusted to concentrate grid points in regions of high curvature. This gives the user consider- 
able flexibility in trying to develop an acceptable grid. An example of the grid generator output is 
shown in Figure 1, the grid for an iced NACA0012 airfoil. 

The equations solved in this grid generator are hyperbolic and take the form; 


x ^ + y^ = 0 e * 1 

"-Vs = 7 = v Eq2 

where x and y are the Cartesian coordinates and £ and q are the coordinates of the body-fitted sys- 
tem. The inverse of the Jacobian, J, approximates the local cell volume, V. Since this is a set of hy- 
perbolic equations, they are solved using a marching procedure starting at the body surface and 

moving along the r\ direction out to the outer boundary. 
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One reason for choosing this grid generator is that it has been modified to work for surfaces with 
concave curvature or surface slope discontinuities. This is especially helpful when trying to re- 
grid the geometry after the ice shape has been added. This grid generator has been used on several 
complex shapes as documented by Barth 7 . It has also been used extensively for ice accretion ge- 
ometries and has produced usable grids, for most cases with minimal user interaction. The issue 
of grid spacing in conjunction with the accuracy of the flow solution can be a problem even for 
clean airfoil geometries, however careful monitoring of the flow solution should allow a high de- 
gree of confidence in the solution. 

Since the overall LEWICE/E code is modular, an alternate grid generator can be easily substituted 
for the current module. The only criterion for substitution is that the grid generator output the grid 
coordinates in PLOT3D format (a NASA developed 3D plotting package). This allows the grid 
file to be used by the other modules as well as by the PLOT3D visualization package. 

n-2. Euler Flow Solution Module 

The Euler code used in LEWICE/E is the inviscid form of the ARC2D code developed by Steger^ 
and Pulliam^. This code uses the grid file and an input file which describes the flow conditions 
and selects some code procedure parameters. The code solves the Euler equations in the body-fit- 
ted coordinate system. The equations solved are. 


d T <2+a^+a/ = 0 


Eq. 3 
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being the contravariant velocities. The Cartesian coordinates (x,y,t) and the body-fitted coordi- 
nates (£,T|,t) are obtained from the grid generator. Details of the numerical algorithm used to 
solve these equations are found in references 8 and 9. 

The flow solver produces velocities and pressures at grid locations on and off the body surface. 
These values are used for two purposes in the LEWICE/E code. They contribute to the evaluation 
of the particle trajectories and to the determination of boundary layer edge conditions used in the 
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integral boundary layer calculation of the energy balance/ice growth module. The user determines 
the validity of the flow field solution by examining residual and lift histories and the pressure co- 
efficient distribution and entropy contours. The flow field solution values are contained in the £ 
matrix, which is stored in PLOT3D format for plotting and to insure a standard file transfer pro- 
cess. 

n-3. Particle Trajectory Module 

Once the flow field for a given geometry has been obtained, the velocities are used to determine 
the trajectories of super-cooled water droplets from release points upstream of the body to impact 
locations on the surface. In the original LEWICE code, this was done with a Lagrangian predic- 
tor-corrector method in a code developed for LEWICE 10 . This code is also used for LEWICE/E, 
however, some changes were required in order to deal with a grid-based flow code instead of the 
original inviscid panel method, which can determine velocity values at any arbitrary location. 

The equations of motion for water droplets are based on the assumption that they are rigid spheres 
and hence the only forces acting are drag and gravity. This assumption is considered valid for 
droplet radii of less than 500 pm 1 . The governing equations are thus, 


m'x = -Dcosy+mgsina 
my = -Dsiny-mgcosa 


Eq. 6 


where m is the droplet mass, a is the angle of attack for the airfoil, and 


y = tan - 


V -V 
l y p y 


x -V 

*p V X 


p y 

D = 


= jd p -v x ) 2 + 0>,-v) 


Eq. 7 

Eq. 8 
Eq. 9 


and where A p is the characteristic area of the droplet, is the drag coefficient of the droplet, p a is 
the density of air, V x and Vy are the components of the flow field velocity at the droplet location 
and x p , y p , x p and y p are the components of the droplet velocity and acceleration, respectively. 

The water droplet trajectories are calculated using an Adams-Boulter predictor-corrector algo- 
rithm. After each calculation of a droplet’s location it is checked to determine if it has impacted on 
the body. The calculation continues until a droplet impacts the surface or until it has passed some 
user selected location, indicating that the droplet has missed the body. Details of this calculation 
are found in reference 1. 


5 



As mentioned in the introduction, the original LEWICE code requires a pseudo-surface in order to 
avoid unrealistic velocities near the surface which lead to inaccurate droplet trajectories. The use 
of a grid based code avoids this problem. Any inviscid flow solution does require some special 
treatment near the surface. The slip flow at the surface is not physically correct and thus the drop- 
let trajectory calculation does not use those values. Instead, once the particle has crossed the first 
grid line off of the surface, the droplet path is no longer altered by the flow field and is considered 
to be tangent to the previously calculated path line. The intersection of that path line with the body 
surface is taken as the impingement location. 

The pattern of droplet impingement on the body surface determines the amount of water that hits 
the surface and becomes part of the ice growth process. The ratio of the actual mass that impinges 
on the surface to the maximum value that would occur if the droplets followed straight-line trajec- 
tories, is called the total collection efficiency, E m . The total collection efficiency for the body is 
found by integrating the local collection efficiency, [J, between the upper and lower limits of drop- 
let impingement. The local collection efficiency is defined as the ratio, for a given mass of water, 
of the area of Impingement to the area through which the water passes at some distance upstream 
of the airfoil. Taking a unit width as one dimension of both area terms, the local collection effi- 
ciency can then be defined as, 



A y 0 

As 


Eq. 10 


where Ay 0 is the spacing between water droplets at the release plane and As is the distance along 
the body surface between the impact locations of the same two droplets. The local collection effi- 
ciency is illustrated in Figure 2. 

The local collection efficiency is the necessary input for the energy balance/ice growth module. It, 
along with the free stream velocity and the cloud liquid water content, determines how much wa- 
ter impinges on the local region of the surface under consideration. Variations in the local collec- 
tion efficiency can significantly alter the ice growth for that surface region. 

n-4. Energy Balance/Ice Growth Module 

The growth of ice on the surface is a complex fluid dynamics, heat transfer, and mass transfer pro- 
cess. The incoming water may freeze on impact or some fraction may freeze while the remaining 
water either runs along the surface or collects in pools. The processes determining which of these 
occur include surface tension effects, roughness, skin friction between water and ice or water and 
airfoil surface, shear forces between water and air, and convection and conduction heat transfer. A 
simplified model of this process has been developed by Messinger 1 1 and is used in the original 
LEWICE code. This model is used in the LEWICE/E model as well. However, as alternate ice 
growth models are developed, they may easily be substituted into this code due to its modular na- 
ture. 

The Messinger model, as implemented in the LEWICE code, is described fully in Appendix A of 
the LEWICE User’s Manual 1 . As such, only a brief description will be given here. The ice growth 
process is modeled as a control volume analysis. The control volume is bounded by the body sur- 
face and by an arbitrary boundary considered to be at the edge of the boundary layer. The two 
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chordwise boundaries coincide with constant % grid lines established by the grid generation code. 
The lower boundary of the control volume is initially the body surface but it moves outward with 
the ice growth process, remaining at the top of the solid-fluid interface. For dimensional com- 
pleteness, the control volume is considered to extend one unit length in the spanwise direction. 

The control volume is used to perform a mass and energy balance, as depicted in Figures 3 and 4. 
The equations governing the mass and energy balance are. 


m. + m. -m.-m. = m 

L r in c r out 


Eq. 11 


where m c is mass flow rate of incoming water, m r is the mass flow rate of runback water from 
the previous control volume, m e is the mass flow rate of evaporated water, m * s mass A° w 
rate of water running back to the next control volume, and is the mass flow rate of water leav- 
ing the control volume due to freeze-out 


and 


T sur (* — 1) 


= ('Mv, sur + *rjw, sur + sur + ^c As + * lA *) 


Eq. 12 


where i w T is the stagnation enthalpy of the incoming water droplets, i w sur ^ _ jj is the enthalpy 
of the water flowing into the control volume from upstream, i v> sur is the enthalpy of the vapor 
leaving the control volume due to evaporation, i w sur is the enthalpy of the water running back to 
the next control volume, i i sur is the enthalpy of the ice leaving the control volume, q c is the heat 
transfer due to convection, and q^ is the heat transfer due to conduction at the bottom of the con- 
trol volume. 

The incoming energy due to water droplet impingement and runback are calculated from known 
information. The energy leaving the control volume due to evaporation and convection can be cal- 
culated independently. The heat transfer due to conduction is not considered in this analysis, as 
the ice layer is considered to act as an insulating surface. This leaves the energy loss due to freeze- 
out and the energy leaving the control volume due to runback as unknowns. In particular, the mass 
flow rates for these two terms are unknown, as was the case for the mass balance. This leaves two 
equations and two unknowns and the system can be solved. The details of the evaluation of each^ 
of the terms in the energy equation can be found in Appendix A of the LEWICE User’s Manual. 

A useful concept for evaluation of the nature of the ice accretion being calculated is the freezing 
fraction. This is the fraction of the total water coming into the control volume that changes phase 
to ice. The equation defining freezing fraction is, 


f= ZL 


m. 


m c + m 


Eq. 13 
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This term can also be used to simplify the evaluation of the energy balance. 

Expanding the terms in the energy equation as described in the LEWICE manual and combining 
Eqs. 12 and 13 yields the following form of the energy equation. 
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+ rh 'J c P. M (7 W 1) " 273 - 15 )l + 


-m e [c pw JT sur -n3.l5) + L v] 
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sur e 
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As 


Eq. 14 


The convection heat transfer term plays an important role in the LEWICE/E energy" balance. It is 
through this term that the aerodynamics and the roughness levels can influence the development 
of the ice accretion. Currently, the convection heat transfer is determined from an evaluation of 
the boundary layer growth on the surface, using an integral boundary layer method. The pressure 
distribution determined by the inviscid panel method is used as input to the boundary layer calcu- 
lation. The boundary layer calculation determines the displacement thickness and the momentum 
thickness. The Reynold’s analogy is used to determine the heat transfer coefficient. Roughness is 
accounted for by a correlation developed by Ruff. 1 The complete description of the integral 
boundary layer calculation is found in Appendix B of the LEWICE User’s Manual. 1 

The solution process for this equation is started by identifying the stagnation point. Since no run- 
back water can enter the control volumes on either side of the stagnation point, one of the un- 
knowns is determined for these two control volumes. Thus, the solution may be marched back, on 
the upper and lower surfaces, towards the trailing edge from the stagnation point. As rh i is deter- 
mined For each control volume, the resulting ice growth thickness can be found from the ice den- 
sity and the dimensions of the control volume. The ice thickness values define a new iced airfoil 
geometry by adding that thickness to the body in a direction normal to the surface. In regions of 
high curvature, the new ice surfaces can intersect or diverge. The method for dealing with the re- 
definition of the surface under these circumstances is the same as that of the original LEWICE 
code and is described in Reference 1. Once the new geometry has been defined, the entire process 
can be started again at the grid generation step. 
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in. Example Cases 

The LEWICE/E code can be used for cases at high subsonic Mach numbers. However, it must 
also be tested at conditions for which the original LEWICE code produced acceptable results. 
Since this effort is not an exhaustive study of the code’s capabilities, but rather an introduction of 
a new approach, two example cases will be described. The first is a reproduction of the LEWICE 
Example 1 calculation for a NACA 0012 airfoil. The second is a calculation for a NACA 65A004 
airfoil at 0° angle of attack. This airfoil was chosen because it was designed for high subsonic 
Mach numbers and because its impingement limits were determined in 1955 by Brun and Vogt 12 , 
using a differential analyzer. 

ffl-1. NACA 0012 

Example case 1 from the LEWICE User’s Manual is a NACA 0012 airfoil at 0° angle of attack 
with glaze icing conditions prevailing. The accretion time is 2 minutes, calculated in two one- 
minute intervals. After the first one minute time interval, the flow field and water droplet trajecto- 
ries are recalculated for the new airfoil geometry, which includes ice on the leading edge of the 
airfoil. The icing conditions are indicated in Table 1. 

The results for the first time step of this case are presented in Figures 5-9, which show the com- 
parison of selected results from the LEWICE and LEWICE/E calculations. Figure 5 shows the lo- 
cal collection efficiency values, plotted as a function of surface distance. The two results compare 
well for impingement limits, however, the LEWICE values at the stagnation region are lower than 
the LEWICE/E values. This is due to the use of the pseudo-surface in the LEWICE calculation. 
The resulting surface is more blunt than the actual surface and results in lower collection efficien- 
cy values. The LEWICE/E calculation uses the actual surface and hence there is less deviation of 
the particle trajectory paths. 

Figures 6 and 7 show the results for the convective heat transfer coefficient. In general, the values 
are the same for both calculations with some differences at the stagnation region. Figure 7 shows 
a close-up of the heat transfer values in that region. Both codes show peaks in heat transfer at 
s=+0.01m from the leading edge. The LEWICE results indicate a drop in the region near the stag- 
nation point, while LEWICE/E shows a decrease with a secondary peak centered on the stagna- 
tion point. This result is due to the differences in velocity values calculated by the two codes in the 
region near the stagnation point. The LEWICE code predicts an edge velocity in this region (i.e. 
±0.002m of the stagnation point) of 43 m/s while the LEWICE/E code predicts a value of 56 m/s. 
These differences in velocity values diminish as the distance from the stagnation region increases. 

Figures 8 and 9 show the results of the ice growth calculation for the first one minute time step. 
The LEWICE/E results show a larger amount of ice near the stagnation region. This is a result of 
the increased amount of water predicted by this analysis coupled with the higher heat transfer in 
that region. The extent of icing along the chord is approximately the same for both calculations, 
reflecting the similarity in collection efficiency away from the stagnation region. 

The second time step continues the calculation for an additional 60 seconds. The new ice shape 
geometry is used to produce a grid, as shown in Figure 1. This grid requires a large number of grid 
points to resolve the leading edge region. The flow field results indicate large pressure spikes at 
the comers of the ice shape. These results are shown in the Cp distributions of Figure 10 and in the 
pressure contour plots of Figure 1 1 . The ability of the code to find a converged solution, using the 
supplied grid, is shown in Figures 12 and 13. Figure 12 shows the residual history for the calcula- 
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tion and indicates that the results converge to machine zero in approximately 1200 iterations. Fig- 
ure 13 shows the entropy contours for the region near the leading edge ice shape. These results 
show that the regions of non-isentropic behavior, indicated by the non-zero contour levels, are 
confined to small spots near sections of high surface curvature and do not seem to contaminate the 
flow field at any distance away from the surface. 

Figures 14-17 show comparisons of the LEWICE and LEWICE/E results for the second time step. 
The collection efficiency values, shown in Figure 14, indicate that the two calculations are quite 
different. This is due to differences in shape resulting from the previous time step. These differ- 
ences, along with the heat transfer coefficient differences shown in Figures 15, combine to pro- 
duce different ice shape results for the two calculations. The resulting ice shapes, shown in 
Figures 16 and 17, differ most markedly at the stagnation region. The LEWICE results produce a 
concave region surrounded by two ice horns. The LEWICE/E results produce a convex region 
around the stagnation point with two less pronounced horns at the edges. Comparison with the ex- 
perimental results, shown in Figure 18, indicate that the LEWICE/E results captured the central 
bulge. The LEWICE results appear to capture the horns somewhat better. Altering the roughness 
parameter could result in lower freezing fractions near the hom regions and thus warrants further 
investigation. 

In general, the results of the two calculations agree with the experimental results well within the 
variability of the ice shape along the span of the airfoil. The LEWICE/E results appear to reflect 
the differences that can result from the elimination of the pseudo-surface. A more exhaustive 
comparison should point out any major differences between the two codes. 

m-2. NACA 65A004 

This case was chosen because the airfoil has a thinner profile, more indicative of high Mach num- 
ber conditions, and because there has been some prior analysis of droplet trajectories for this air- 
foil. The conditions for this calculation are listed in Table 2. The relatively high Mach number of 
this case (i.e. M = 0.58) was used in order to highlight the capabilities of a compressible flow 
analysis. 

In this case, the LEWICE/E calculations were compared to an earlier analysis performed by Brun 
and Vogt 12 . They calculated collection efficiencies for this airfoil at several angles of attack and 
several combinations of velocity and particle size. Only one case was examined for this investiga- 
tion. However, a more exhaustive analysis would be useful and is contemplated for a future study. 

The collection efficiency results from the two analyses were compared and are shown in Figure 
19. The two methods agree very well on the overall shape of the curve. The impingement limits 
differ by approximately 0.01 meters. This difference amounts to about 0.5 percent of the chord 
length. Brun and Vogt indicate that their results near the stagnation region are not as accurate as 
further aft along the airfoil surface. In any case, the LEWICE/E results indicate a peak (3 value of 
0.89 at the stagnation point, while the peak [3 value of Brun and Vogt is hard to estimate from the 
graph presented in their paper. 

The LEWICE/E code was also used to investigate a range of temperature conditions (i.e. 267-242 
K) for this airfoil while keeping the other icing conditions constant. The results are presented in 
Figures 20 and 21, showing the ice profile and freezing fraction results respectively. The ice pro- 
file results show the normal progression from running wet through glaze to rime conditions. The 
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freezing fraction results also reflect this progression, going from zero everywhere at the warmest 
condition to one everywhere at the coldest condition. 

There is an interesting result for a warm temperature condition, shown in Figure 20b. For this cal- 
culation, the static temperature was -1 1°C. The ice shape profile indicates no growth near the 
stagnation point followed by a region of ice growth further aft along the surface. This behavior is 
also seen in the freezing fraction results, shown in Figure 21b. The freezing fraction near the stag- 
nation region is zero and then at some point along the surface it increases until beyond the im- 
pingement limits. 

The free stream stagnation temperature for this case is 7°C. Thus, over the first portion of the sur- 
face the convective heat transfer is heating the surface and preventing freezing. Apparently, the 
stagnation temperature variation along the surface is such that there is some point at which the 
evaporative cooling balances the convective heating. Once this occurs the water must start to 
freeze. This condition deserves more attention and will be studied further in a later investigation. 

The roughness value used for this case was 0.00286 meter. This value is consistent with the guide- 
lines suggested in the LEWICE manual. However, when compared to the leading edge radius for 
this airfoil, which is 0.002 meter, it appears to be inordinately large. Work by other authors ’ has 
suggested that the roughness parameter may require some modification. In the case of thin air- 
foils, such as the NACA 65A004, the roughness correlation of LEWICE may also require further 
study. 

IV. Conclusion 

An Euler flow analysis method has been combined with the trajectory and ice growth prediction 
calculation routines of the LEWICE code to produce a compressible flow ice accretion code, 
LEWICE/E. The LEWICE/E code was compared to the original LEWICE code for example case 
1 from the LEWICE User’s Manual. Results indicate general agreement between the two codes 
and reasonable reproduction of the experimental results. Differences between the codes did 
emerge in (5 values and convective heat transfer coefficients. The collection efficiency differences 
are due to the fact that LEWICE/E does not use a pseudo-surface in the trajectory calculation. The 
heat transfer differences are a result of the differing velocity values obtained for the stagnation re- 
gion by the two flow field codes. The LEWICE/E code was also used to examine a thin airfoil pro- 
file at a Mach number of 0.58. Temperature sweep results indicated the normal progression from 
glaze to rime conditions as the temperature decreased. 

The results of this development activity suggest the course of future work. A parameter study 
needs to be undertaken to determine which icing conditions are most sensitive to compressibility 
effects. The role of compressibility in the development of the integral boundary layer equations 
requires further study, especially in regard to the convective heat transfer coefficient. The sensitiv- 
ity of ice shape prediction to grid resolution should also be investigated. 

Finally, the Euler code used in this study can also be run as a Navier-Stokes code. Thus, the next 
step in the development of the LEWICE code will be to incorporate Navier-Stokes calculation 
procedures. The Navier-Stokes code will be used to replace the integral boundary layer calcula- 
tion in the ice growth routines and to produce lift and drag values for iced airfoils. This will result 
in a complete icing analysis capability for two-dimensional geometries. 
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Example 1 : NACA 0012 Airfoil 

Velocity 

= 129.46 m/s 

Static Temperature 

= 260.55 K 

Static Pressure 

= 90748.0 Pa 

LWC 

= 0.50 g/m 3 

Droplet Diameter 

= 20.0 jim 

Roughness Height 

= 0.00035 m 

Chord Length 

= 0.30 m 


TABLE 1 : ICING CONDITIONS FOR EXAMPLE CASE 1 


Example 2 : NACA 65A004 Airfoil 

Velocity 

= 190.13 m/s 

Static Temperature 

= 242 - 267 K 

Static Pressure 

= 90900.0 Pa 

LWC 

= 1.50 g/m 3 

Droplet Diameter 

= 25.0 |im 

Roughness Height 

= 0.00286 m 

Chord Length 

= 1.96 m 


TABLE 2 : ICING CONDITIONS FOR EXAMPLE CASE 2 
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FIGURE 1 - NEAR FIELD GRID FOR NACA0012 AIRFOIL WITH 

ICE SHAPE ON LEADING EDGE. 



FIGURE 2 - DEFINITION OF LOCAL COLLECTION EFFICIENCY. 
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FIGURE 3 - MASS BALANCE CONTROL FIGURE 4 - ENERGY BALANCE CONTROL 

VOLUME. VOLUME. 



S (U) 


FIGURE 5 - BETA DISTRIBUTIONS FOR NACA 0012 AIRFOIL. 
COMPARISON OF LEWICE AND LEWICE/E CALCULATIONS. 
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HTC 



FIGURE 6 - HEATTRANSFER COEFFICIENT DISTRIBUTIONS FOR NACA 
0012 AIRFOIL. COMPARISON OF LEWICE AND LEWICE/E 



S (m) 


FIGURE 7 - HEAT TRANSFER COEFFICIENT DISTRIBUTIONS NEAR 
LEADING EDGE FOR NACA 00 12 AIRFOIL. COMPARISON OF LEWICE 
AND LEWICE/E CALCULATIONS. 
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FIGURE 9 - LEWICE/E ICE SHAPE PREDICTION FOR NACA 0012 AIRFOIL 
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FIGURE 10- PRESSURE COEFFICIENT FIGURE 11- PRESSURE CONTOURS NEAR 

DISTRIBUTION ON ICED NACA 0012 AIRFOIL. THE LEADING EDGE OF AN ICED NACA 

0012 AIRFOIL. 



FIGURE 12 - RESIDUAL HISTORY FROM 
LEWICE/E CALCULATION FOR ICED NACA 
0012 AIRFOIL 



FIGURE 13 - ENTROPY CONTOURS NEAR 
THE LEADING EDGE OF AN ICED NACA 
0012 AIRFOIL. 
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I.Or 



FIGURE 14 - BETA DISTRIBUTIONS FOR ICED NACA 0012 AIRFOIL. 



FIGURE 15 - HEAT TRANSFER COEFFICIENT DISTRIBUTIONS NEAR LEADING EDGE FOR ICED NACA 0012 
AIRFOIL COMPARISON OF LEWICE AND LEWICE/E CALCULATIONS. 
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FIGURE 1 8 - EXPERIMENTAL ICE SHAPE TRACING FROM NACA 0012 
AIRFOIL. TIME = 120 SEC. (REF. 1) 



FIGURE 19 - BETA DISTRIBUTION FOR NACA 65A004 AIRFOIL. COMPARISON OF LEWICE/E 
CALCULATIONS AND ANALYTICAL RESULTS OF BRUN & VOGT. 
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